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Abstract 

We introduce Gaussian Matrix Product States (GMPS), a general- 
ization of Matrix Product States (MPS) to lattices of harmonic oscil- 
lators. Our definition resembles the interpretation of MPS in terms 
of projected maximally entangled pairs, starting from which we de- 
rive several properties of GMPS, often in close analogy to the finite 
dimensional case: We show how to approximate arbitrary Gaussian 
states by MPS, we discuss how the entanglement in the bonds can 
be bounded, we demonstrate how the correlation functions can be 
computed from the GMPS representation, and that they decay ex- 
ponentially in one dimension, and finally relate GMPS and ground 
states of local Hamiltonians. 

[This work originally appeared as Sec. VII of quant-ph/0509166, and 
is published in the Proceedings on the conference on Quantum infor- 
mation and many body quantum systems, edited by M. Ericsson and 
S. Montangero, pg. 129 (Edizioni della Normale, Pisa, 2008).] 



1 Introduction 



In the last years, there has been considerable activity on the border between 
quantum information theory and condensed matter physics, and quantum 
information concepts have successfully been applied to the description of 
quantum many-body systems. An important step has been the interpre- 
tation of Matrix Product States (MPS) in terms of projected maximally 
entangled pairs. MPS form a hierarchy of states which prove very success- 
ful as a variational ansatz for simulating ground states of one-dimensional 
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quantum systems, as done in the Density Matrix Renormalization Group 
(DMRG) method [HI2]- From the perspective of quantum information, MPS 
are formed by taking virtual maximally entangled pairs between adjacent 
sites and applying a linear map on each site to obtain the physical system [3] . 
This entanglement-based description led to a better understanding of MPS 
and gave rise to new algorithms and extensions of DMRG to e.g. thermal 
states, time evolutions, and higher dimensional systems [IH7], but also to 
new analytical tools for investigating e.g. quantum phase transitions, renor- 
malization group transformations, or the sequential generation of quantum 
states [MO]. 

Given the success of the MPS framework in the description of finite- 
dimensional spin systems, it is natural to look for generalizations to e.g. 
bosonic or fermionic systems. In this paper, we introduce bosonic Gaussian 
Matrix Product States (GMPS), which describe Gaussian states on lattices 
of harmonic oscillators (i.e., bosonic modes). Such systems are frequently 
realized in physical setups, e.g. by the vibrational modes of ions in linear 
traps or by arrays of nanomechanical oscillators, and since they are typically 
goverened by quadratic Hamiltonians, their ground and thermal states are 
Gaussian. 

Our definition of GMPS resembles the quantum information perspective 
on MPS, where one takes maximally entangled pairs and applies a linear map 
to obtain the physical system. Starting from this definition, we show that ev- 
ery (translation invariant) Gaussian state can be represented as a (translation 
invariant) GMPS, and discuss how to minimize the amount of entanglement 
used in the bonds - different from the finite-dimensional case, this is an issue 
since bosonic bonds can carry an unbonded amount of entanglement. We 
discuss the properties of two-point correlation functions of GMPS and show 
that they can be easily computed from the GMPS representation; for the 
case of pure one-dimensional GMPS, we prove that the correlations decay 
exponentially (as it is the case in finite dimensions) and explicitly derive the 
correlation length. We end our discussion on GMPS by showing that - again 
in analogy to the finite dimensional case - every GMPS is the ground state 
of a local Hamiltonian. 

Since Gaussian states are completely characerized by their second mo- 
ments and thus by a number of parameters quadratic in the system size, 
unlike for spin systems Gaussian MPS will not by themselves yield an ex- 
ponentially more efficient parametrization. However, they can be used to 
describe translational invariant states with a constant number of parame- 
ters and thus also in the limit of an infinite chain. Moreover, the GMPS 
parametrization should have favorable properties e.g. for variational mini- 
mizations with respect to local observables. 
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The results presented in this work have already been applied in [TTj[T2]. 



2 Gaussian states 

Consider a system of N bosonic modes which are characterized by N pairs 
of canonical operators (Q±, Pi, . . . ,Qn, Pn) ='■ R, and where the canonical 
commutation relations (CCR) are governed by the symplectic matrix a via 

JV 



( 1 \ 
[R k ,R t ] = ia M , cr = 0( 1 Q J 

n=l ^ ' 



Then, Gaussian states are defined as states which have a Gaussian Wigner 
distribution in phase space. Those state are frequently met in physics, since 
ground or thermal states of quadratic Hamiltonians are Gaussian states, and 
evolution under a quadratic Hamiltonian leaves Gaussian states Gaussian. 
They are completely characterized by their first moments d k = tr[pRk\ 
(which can be changed by local operations, so that we set them to zero 
w.l.o.g.) and their covariance matrix (CM) 



Iki = tr 



p{R k - d k ,Ri - di] 



X 



where {-, •}+ is the anticommutator. The CM satisfies 7 > ia, which ex- 
presses Heisenberg's uncertainty relation and is equivalent to the positivity 
of the corresponding density operator p > 0. Purity of the state is charac- 
terized by det7 = 1 or equivalently (cy) 2 = — 1. 

When the state under consideration is translational invariant, it often 
proves convenient to describe it in the Fourier basis. For simplicity, let us 
consider a one- dimensional translational invariant chain of length N with 
periodic boundaries and reorder the canonical operators such that R = 
(Qt, Q N , Pi, ... , P N )- Then, 

( 1Q 1QP 
\ Tqp 1p 

and translational invariance is reflected by the fact that any matrix element 
A k i (where A G {7q,7p,7qp}) depends only on the distance k — I (which is 
understood mod N), so we can write A k _i := A k i. Matrices of this type are 
called circulant and are simultaneously diagonalized by the Fourier trans- 
form. We write the diagonal elements of the Fourier transform as a function 
of the angle <fi = 2mn/N for m = 0, . . . , N — 1; the Fourier transform of a 
cirulant matrix A then reads A(4>) = ^A n e~ m ^. 
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An interesting property of translational invariant pure states with one 
mode per site is that they are point symmetric JT3] . This can readily be seen 
from the representation 

^ = ( 1Q 1QP \ = ( X XY \ 

\1qp IP J \ YX X~ l + YXY ) 

in Q-P partitioning with X > 0, Y = Y T real [13], since X and Y have to 
be circulant and therefore commute. Hence, jqp = XY = YX = 7q P , i.e., 
7 is point symmetric. Note that this implies in particular that the Fourier 
transform 7(0) is real. 

3 Gaussian Matrix Product States 

In the following, we introduce Gaussian matrix product states (GMPS). The 
definition resembles the physical interpretation of finite-dimensional matrix 
product states as projected entangled pairs: In finite dimensions, MPS can be 
described by taking maximally entangled pairs of dimension D between adja- 
cent sites and applying arbitrary local operations on each site, thus mapping 
the D x D dimensional input (the virtual system) to a <i-dimensional output 
state (the physical system). Similarly, GMPS are obtained by taking a num- 
ber of entangled bonds and applying local (not necessarily trace-preserving) 
operations T^, where the boundary conditions can be taken either open or 
closed. Any GMPS is completely described by the type of the bonds and 
by the operations T^. Note that this construction holds independent of the 
spatial dimension. For one dimension, it is illustrated in Fig. [TJ As ma- 
trix product states are frequently used to describe translationally invariant 
systems, an inportant case is given if all maps are identical, = T Vi. 



i : i+1 




Figure 1: Construction of Gaussian Matrix Product States (GMPS). GMPS 
are obtained by taking a fixed number M of maximally entangled (i.e., EPR) 
pairs shared by adjacent sites, and applying an arbitrary 2M to 1 mode 
Gaussian operation on site i. 
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Figure 2: The Jamiolkowski isomorphism. The Gaussian channel described 
by the state T can be implemented by projecting the input state 7; n (mode 
A) and the input port of T (mode B) onto the EPR state (symbolized by 
curly brackets). In case of success, the output is obtained in mode C. The 
operation can be made trace-preserving by measuring in a basis of displaced 
EPR states, and displacing C according to the measurement outcome. 

In order to define MPS in the Gaussian world, we have to decide on the 
type of the bonds as well as on the type of operations. We choose both the 
bonds to be Gaussian states and the operations to be Gaussian operations, 
i.e., operations mapping Gaussian inputs to Gaussian outputs. For now, we 
will take the bonds to be maximally entangled (i.e., EPR) states, such that 
the only parameter originating from the bonds is the number M of EPRs. 
We show later on how the case of finitely entangled bonds can be easily 
embedded. 

As to the operations, we will allow for arbitrary Gaussian operations. 
Operations of this type are most easily described by the Jamiolkowski iso- 
morphism [15]. There, any Gaussian operation T which maps N input modes 
to M output modes can be described by an N + M mode covariance matrix V 
with block B (input) and C (output). The corresponding map on some input 
state 7; n in mode A is implemented by projecting the modes A and B onto an 
EPR state as shown in Fig. [21 such that the output state T(jm) is obtained in 
mode C. Conversely, the matrix T which represents the channel T is obtained 
by applying the channel to one half of a maximally entangled state. The du- 
ality between T and V is most easily understood in terms of teleportation, 
and shows that this characterization encompasses all Gaussian operations. 
Note that the protocol of Fig. [2] can be always made trace-preserving by pro- 
jecting onto the set of phase-space displaced EPR states and correcting the 
displacement of mode C according to the measurement outcome [16] . 

In the following, we will denote all maps T by their corresponding CM V. 
Sometimes, we will speak of the modes B and C as input and output ports 
of T, respectively. 

We now discuss how the covariance matrix of the output will depend on 
the CM of the input and on the channel V [16l[T7] . This is most easily com- 
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puted in the framework of characteristic functions [18]. The characteristic 
function of the output is given by 

Xc{&) oc J e-^r^A e -tl c rt B c 5 ( XA _ XB } 5 ( pA+pB ^ AB ? 
and by integrating over subsystem A, we obtain 

Xc{&) oc J e-^c M ^c^ B 

with 



M 



Basically, the integration J (1£a5(xa — xb)8(pa+Pb) does the following: first, 
it applies the partial transposition 9 = ( J _?i ) to one of the subsystems, and 
second, it collapses the two systems A and B in the covariance matrix by 
adding the corresponding entries. The integration over on the other 
hand, leads to a state whose CM is the Schur complement of Mix, M22 — 
M 2 iM^ 1 1 M12, such that the output state is described by the CM 

7out = Tc — ^cbt; — — 7. — tFbC ■ 
+ 6>7 in 6> 

Let us briefly summarize how to perform projective measurements onto 
the EPR state in the framework of CMs, where we denote the measured 
modes by A and B, while C is the remaining part of the system. First, apply 
the partial transposition to B, second, collapse A and B, and third, take the 
Schur complement of the collapsed mode AB, which gives the output CM of 
C. 

In analogy to the finite-dimensional case, we will focus on pure GMPS. 
Particularly, a GMPS is pure if the which describe the operations are 
taken to be pure, which we assume from now on. Let us finally emphasize 
that the given defintion of MPS holds independent of the spatial dimension 
of the system, as do most of the following results, and in fact applies to an 
arbitrary graph. 



4 Completeness of Gaussian MPS 

In the following, we show that any pure and translational invariant state can 
be approximated arbitrarily well by translational invariant Gaussian ma- 
trix product states, i.e., GMPS with identical local operations T. (Without 
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translational invariance, this is clear anyway: the complete state is prepared 
locally and teleported to its destination using the bonds.) The proof is pre- 
sented for one dimension, but can be extended to higher spatial dimensions. 
Note that a similar result also holds for finite-dimensional MPS |19j . 

Given a translational invariant state 7, there is a translational invariant 
Hamiltonian H which transforms the separable state 1 into 7, 7 = SS T , 
S = e aH . It has been shown [2D] that this time evolution can be approximated 
arbitrarily well by a sequence of translational invariant local (one-mode) and 
nearest neighbor (two-mode) Hamiltonians Hj, 

j 

e^«IJe©»^, (2) 

j=i 

where the Hj act on one or two modes, respectively, and approach the identity 
for growing J. 

Clearly, translational invariant local Hamiltionians can be implemented 
by local maps without using any EPR bonds. In the following, we show how 
translational invariant nearest-neighbor interactions can be implemented by 
exploiting the entanglement of the bonds. The whole procedure is illustrated 
in Fig. [3] and requires two EPR pairs per site. We start with some initial 
state 7; n onto which we want to apply = e® aHj ~ 1 + @ n crHj. 

First, we perform local EPR measurements between the modes of 7; n and 
one of the bonds in order to teleport the modes of 7i n to the left, cf. Fig. [3k. 
Then, the infinitesimal symplectic operation S = e aU] is applied to the left- 
teleported mode and the second bond, Fig. [3b. In the last step, another 
EPR measurement is performed which teleports the left-teleported mode 
back to the right, and "into" the mode on which the adjacent S was applied. 
As the operations e aHi > ~ 1 + oHj all commute, the "nested" application 
of the nearest neighbor symplectic operations S indeed give S 1 ®, and thus 
the remaining mode really contains the output 7 out = S^inS^. The whole 
decomposition (T2J) can be implemented by iterated application of the whole 
protocol of Fig. 

5 GMPS with finitely entangled bonds 

Let us now consider the entanglement contained in the bonds and show 
that infinitely entangled bonds can be replaced by finitely entangled ones. 
Intuitively, this should be possible whenever the channel destroys some 
of the entanglement of the bond anyway, i.e., is non-maximally entangled. 
In that case, it should be possible to use a less entangled bond while choosing 
a channel which does not destroy entanglement any more. 
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Figure 3: Implementation of 
a translational invariant nearest 
neighbor Hamiltonian in a trans- 
lational invariant fashion. Start- 
ing from 7i n , the input is first 
teleported to the left, then, the 
infinitesimal time evolution S = 
e uH , H <C 1, is performed, and fi- 
nally, the state is teleported back. 



The method is illustrated in Fig. HI Again, for reasons of clarity we restrict 
to one dimension and one bond. The argument however applies independent 
of the spatial dimension and the number of bonds. The only restriction we 
have to make is the restriction to pure GMPS, i.e., those with pure T^. 

Consider a GMPS with local channels given by and infinitely entan- 
gled bonds, Fig. UK First, apply a Schmidt decomposition (21] to in 
the partition A\BC, which can be always done as long as is pure. The 
Schmidt decomposition allows us to rewrite the state as shown in Fig. |4b - 
an entangled state between modes A and C with two-mode squeezing r^, 
B in the coherent state 1, and sympectic operations and Sg C which are 
applied to modes A and BC, respectively. As the bond itself is infinitely 
entangled, we can teleport the sympectic operation through the bond to the 
next site as indicated in Fig. Hb. Then, S^~^ can be merged with Sg C to a 
new operation acting on modes B and C of site i (Fig. 3b). Finally, in the 
triples consisting of one maximally entangled state, one non-maximally en- 
tangled state, and the projection onto the EPR state, the maximally and the 
non-maximally entangled state can be swapped, resulting in Fig. There, 
we have finitely entangled bonds, while the infinite entanglement has been 
moved into the new maps T^. 

It is tempting to apply this construction to the completeness proof of 
the preceding section in order to obtain a construction which is less wasting 
with respect to resources. However, for any iterative protocol this is most 
likely difficult to achieve. The reason for this is found in the no-distillation 
theorem which states that with Gaussian operations, it is not possible to 



S 



Figure 4: How to make the bonds of GMPS finitely entangled, a) The 
initial MPS. b) Do a Schmidt decomposition of the original map T. c) Move 
the through the infinitely entangled bond to the next site, d) Swap the 
finitely and the infinitely entangled pair. 

increase the amount of entanglement between two parties [16]. Particularly, 
this implies that in each step of an iterative protocol, the bonds need to have 
at least as much entanglement as can be obtained at the output of this step, 
maximized over all inputs where the entanglement is increased. This is indeed 
a severe restriction, although it does not imply the impossibility of such a 
protocol. One could, e.g., create a highly entangled state in the first step and 
then approach the desired state by decreasing the entanglement in each step. 
Still, it seems most likely that a sequence of MPS which approach a given 
state efficiently will have to involve more and more bonds simultaneously and 
thus cannot be constructed in an iterative manner. 

6 Correlation functions of Gaussian MPS 

In this section, we show how to compute correlation functions from the maps 
rW which describe the GMPS. We show that this can be done efficiently, 
i.e., in a time which is polynomial in the systems size, independent of the 
dimension of the graph. This is different from the finite dimensional case, 
where correlation functions of e.g. two-dimensional MPS cannot be computed 
efficiently (22]. Of course, this is not too surprising given that Gaussian states 
can be fully characterized by a number of paramaters quadratic in the number 
of modes. 

Let us start with the general case of different T^, as in Fig. \5fc,. The 
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Figure 5: If the local operations are described by states via the Jami- 
olkowski isomorphism, the construction of GMPS can be simplified by re- 
placing the measurement-bond-measurement triples by a simple projection 
onto the EPR state. 



calculation can be facilitated by the simple observation that the triples con- 
sisting of two projective measurements and one EPR pair can be replaced 
by a single projection onto the EPR state, Fig. \5jp. It follows that we can 
apply the formalism for projective measurements onto the EPR state which 
we presented in Section [3j Starting from we first partially transpose 

all B modes, then collapse Ai + \ and Bi for all i, and finally take the Schur 
complement of the merged mode. In case of periodic boundary conditions, 
this can be expressed by the transformation matrix 



n 



i A ne B o 
o o i c 



(3) 



which maps ABC onto A'C, where 9b = 9® 1 is the partial transposition on 
system B, and 1Z is the circulant right shift operator, (71) ij = mod jv®!- 

Then, the output state, i.e., the GMPS characterized by rW, is 



7 = SC A / 




where SCx(U) is the Schur complement of the X part of U, SCx(U) = 
Uyy — UyxUxxUxy ■ For open boundary conditions, the matrix II has to 
be modified accordingly at the boundaries. All the involved operations scale 
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polynomially in the product NM of the number of sites N and the number 
of modes M. 

In case all the local maps are chosen equal, = r Vi, the above formula 
can be simplified considerably. Therefore, note that the Fourier transform 
can be taken into the Schur complement, and that II as well as ® i=1 = 
r ® Ijv are blockwise circulant so that both are diagonalized by the Fourier 
transform. Thereby, r ® 1 is mapped onto the constant function T, and the 
same holds for 1 and 9 in ([3]). The right shift operator 1Z, on the other hand, 
is transformed to e %< ^\: the EPR measurement performed between adjacent 
sites leads to a complex phase of (p. Altogether, we have 

nrfr 1 . 

Directly expressed in terms of the map T, this reads 

j(<f>) = T C - r c{AB A T — r At T ABlc (4) 

where A = (1^ ; e l ^6 B ) is the upper left subblock of ft. 



n 



iL A e" v u B u 
In 
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7 States with rational trigonometric functions 
as Fourier transforms 

Let us now restrict to pure MPS (i.e., those for which T is pure) with one 
mode per site. As we have shown in Section [2J those states have reflection 
symmetry and therefore j(<t>) = 7o + 2 J2 n >o ^ n cos ( n 0) is rea l- This implies 
that the sines in (j4j) can only appear in even powers sin 2 ™ = (1 — cos 2 (p) n . 
Therefore, the Fourier transform 7 of any pure Gaussian MPS, which is a 
2x2 matrix valued function of 0, has elements which are rational functons of 
cos(0), {l{4>))xy = Pxy (cos {(j)))/q X y (cos ((f))) with p, q polynomials. The degree 
of the polynomials is limited by the size of ATab^ , and thus by the number 
M of the bonds. One can easily check that dimp < 2M + 1 and dimg < 2M. 

For the following discussion, let us write those rational functions with a 
common denominator d, 

Urh\ - 1 f 9( cos (^)) r(cos(0)) \ . , 

1W} ~ d(cos(0)) V r(cos(0)) p(cos(0)) ) ' ^ 

where q, p, r, and d are polynomials of degree L. Then, the set of all 
such 7 with L > 2M + 1 encompasses the set of translational invariant 
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GMPS with M bonds. Computing correlation functions in a lattice of size N 
can be done straightforwardly in this representation by taking the discrete 
Fourier transform of 7(0) which scales polynomially with N, and in the 
following section we show that for one dimension, the correlations can be 
even computed exactly in the limit of an infinite chain. 

It is interesting to note that 7(0) is already determined up to a finite 
number of possibilities by fixing r and d. Since 7 is pure, 1 = det7 = det'j, 
and therefore, pq = d 2 + r 2 . Therefore, the zeros of pq are the zeros of d 2 + r 2 , 
such that the only freedom is to choose how to distribute the zeros on p and 
q. On the contrary, fixing only q and d does not give sufficient information, 
while choosing p, q and d (i.e., the diagonal of 7) does not ensure that there 
exists a polynomial r such that pq — r 2 = d 2 . 

From the above, it follows that 2L+1 parameters are sufficient to describe 
7(0), where L is still the degree of the polynomials. This encloses all trans- 
lational invariant Gaussian MPS with bond number M < (L — l)/2, which 
need (2M + l)(2M + 2) = L(L + 1) parameters. Therefore, the class of states 
where 7(0) is a rational function of cos(0) is a more efficient description of 
translationally invariant states than Gaussian MPS are. 

Let us stress once more that the results of this section hold for arbitrary 
spatial dimension. 

8 Correlation length 

In the following, we show that the correlations of one-dimensional GMPS 
decay exponentially, and we explicitly derive the correlation length. The 
derivation only makes use of the representation (J5]) of Gaussian MPS and 
thus holds for the whole class of states where the Fourier transform is a 
rational function of the cosine. We will restrict to the case where the state 
T associated to the GMPS map has no diverging entries, which corresponds 
to the case where the denominator d(cos(0)) in ([5]) has no zero on the unit 
circle Q 

The correlations are directly obtained by back-transforming the elements 
of 7(0), which are rational functions [7(0)] s = s(cos(0))/d(cos(0)), s G 
{p, q, r}; in the limit of an infinite chain, 



1 The case where d has zeros on the unit circle corresponds to critical systems, which 
is why the correlations diverge. In the case of a Hamiltonian H = V © 1, however, the 
ground state correlations of P do not diverge [T3]. As in that case one has pq = d 2 , 
p/d= d/q need not have a singularity just because q/d has one. 




12 



Now transform s, d to complex polynomials via cos(</>) — Y (z + l/z)/2, and 
expand with z K , s(z) := z K s(z), d(z) := z K d(z), where K is chosen large 
enough to make s, d polynomials in z. Then, 



s[z)z 



n-l 



E 



d(z) 

1 



-dz 



z;:d(z;)=0 



1)! dz"*- 1 



S(Z)Z 



n-l 



dAz) 



by the calculus of residues, where Vi is the order of the zero Zi in d and 
d~i(z)(z — Zi) Vi = d(z). For n > z/j, Di oc zf~ u '\ and it follows that the 
correlations decay exponentially, where the correlation length is given by the 
largest zero of q(z) inside the unit circle. 

This proof only holds for one-dimensional GMPS. However, it can be 
proven for arbitrary spatial dimensions that the correlations decay faster 
than any polynomial by iterated integration by parts with respect to one 
component of 0, cf. [T3] . 



9 Gaussian MPS as ground states of local 
Hamiltonians 

Finally, let us focus on the relation of translational invariant Gaussian MPS 
and local Hamiltonians. We prove that every GMPS is the ground state 
of a local Hamiltonian, while conversely most Hamiltonians do not have 
GMPS as an exact ground state - again, this is in close analogy to the finite- 
dimensional case P35]. Once more, the proof only requires the state to be 
of the form Eq. ([5]). We will make use of some results on ground states 
of translational invariant quadratic Hamiltonians presented in |13j . Define 
the Hamiltonian matrix H > via the Hamilton operator "H by virtue of 

"H = Yl HkiRkRi, as well as the spectral function £ = \/det H. Then, the 
ground state is given by 

7 = (£@Sy l oHa T (6) 

and has energy \ti£. 

Given a pure state 7 with Fourier transform ()5]), define 

fj(j,\ _ ( P(cosO)) -r(cos(») \ ( . 

W) \ -r(cos(0)) g(cos(0)) J { ] 
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Then, H corresponds to a local Hamiltonian - the interaction range is the 
degree of p,q,r - and £(4>) = [y/pq — r 2 ] (cos <$) = d(cos 0), which together 
with (0) proves that 7 is the ground state of H . 

Let us also have a brief look at the converse question: Given a local 
Hamiltonian, when will it have a GMPS as its ground state? Any local and 
translational invariant Hamiltonian has a Fourier transform which consists 
of polynomials in cos(0), and thus we adapt the notation of Eq. ([7]). Then, 
following Eq. ([6]) the ground state is represented by a rational function of 
cos(0) in Fourier space exactly if pq — r 2 = d 2 is the square of another 
polynomial. In terms of the original Hamiltonian, this implies that HqHp — 
Hq P has to be the square of another banded matrix. For example, for H = 
V © 1 one would need V = X 2 with X again a banded matrix 
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